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ABSTRACT 


In  the  course  of  investigating  the  effects  of  a  nuclear  explosion  on 
electromagnetic  wave  propagation,  it  is  often  convene to  consider 
the  ionized  constituents  of  the  atmosphere  to  be  represented  by  a 
single  species  of  positive  ions,  a  single  negative  ion  species,  and 
electrons.  The  number  densities  of  these  charged  particles  at  a 
specific  time  after  detonation  are  then  obtained  from  the  solution  of 
a  system  of  three  differential  equations  which  describe  the  time 
rate  of  change  of  each  particle  density.  An  exact  analytic  solution 
has  not  been  found  for  these  equations,  and  hence  either  approximate 
analytic  forms  or  numerical  integration  methods  are  employed.  The 
latter  offers  high  accuracy  but  is  often  too  time-consuming  for 
general  application.  In  this  report  a  new  set  of  approximate  analytic 
solutions  are  derived  which  provide  accuracies  approaching  that  of 
numerical  integration  for  many  cases,  but  which  retain  the  advantage 
of  rapid  computation. 
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INTRODUCTION 


An  essential  part  of  an  investigation  of  the  effects  of  nuclear  ex¬ 
plosions  on  electromagnetic  wave  propagation  is  the  determination 
of  atmospheric  ionization  produced  by  energy  released  by  the  deto¬ 
nation.  Once  the  level  of  initial  ionization,  produced  by  the  almost 
instantaneous  energy  pulse,  and  the  continuing  rate  of  ion  production, 
resulting  from  beta  and  gamma  radiation  accompanying  the  radio¬ 
active  decay  of  the  fission  debris,  have  been  determined  from  the 
burst  parameters  and  atmospheric  transmission  properties,  the 
problem  becomes  one  of  obtaining  a  solution  to  the  system  of 
differential  equations  which  describe  the  time  rate  of  change  of  the 
number  density  of  each  atmospheric  constituent.  Well  over  a 
hundred  different  reactions  involving  a  dozen  or  more  different 
atmospheric  species  may  be  included  in  the  general  problem. 

Keneshea  [Ref.  l]  has  investigated  this  problem  and  has  developed 
a  digital  computer  program  which  employs  numerical  integration 
methods  to  treat  a  system  of  thirteen  variable  species  involving 
131  reactions , 

When  many  different  altitudes  must  be  considered  under  varying 
conditions  of  initial  ionization  and  ionization  source  strength,  it  is 
often  not  feasible  to  consider  the  general  problem,  and  some  means 
of  approximating  the  physical  situation  with  a  simpler  mathematical 
model  is  needed.  Furthermore,  many  of  the  coefficients  which 
specify  the  rates  at  which  the  various  reactions  proceed  are  poorly 
known  at  best,  thus  reducing  the  need  for  a  rigorous  treatment  of 
the  problem  in  many  investigations.  One  approach  which  is  widely 
used  consists  of  employing  three  variable  species  to  represent  the 
ones  actually  present:  a  single  positive  ion  species,  a  single 
negative  ion  species,  and  electrons.  The  neutral  population  is  con¬ 
sidered  to  be  constant  at  any  given  altitude.  These  three  species 
are  allowed  to  participate  in  four  reaction*!:  recombination  of  electrons 
with  positive  ions,  mutual  neutralization  o^  positive  and  negative 
ions,  attachment  of  electrons  to  neutrals,  and  detachment  of  electrons 
from  negative  ions.  As  in  the  more  general  approach,  charge  balance 
is  maintained  (the  positive  ion  density  equals  the  sum  of  the  electron 
and  negative  ion  densities). 
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The  problem  is  therefore  reduced  to  a  system  of  three  simultaneous, 
first-order,  nonlinear  differential  equations.  Even  this  reduced 
system  is  still  difficult  to  solve  analytically,  and  resort  must  be 

I  •  jfc 

made  to  either  numerical  methods  or  approximate  solutions.  When 
a  large  number  of  calculations  must  be  made,  approximate  analytic 
solutions  are  often  preferred  for  their  relative  simplicity  and  ease 
of  computation.  The  disadvantage  of  such  solutions  has  been  the 
limited  accuracy  of  the  results  obtained  from  them  under  certain 
combinations  of  parameters.  It  is  an  unfortunate  fact  that  most  of 
the  approximate  solutions  exhibit  their  poorest  accuracy  in  cases 
involving  altitudes  between  roughly  50  krn  and  80  km,  which  is  pre¬ 
cisely  where  the  need  for  accuracy  is  the  greatest  since  most  of  the 
propagation  distunbances  due  to  a  nuclear*  detonation  are  maximum  at 
t  these  altitudes. 

i 

s 

j 

I  The  mathematical  difficulty  which  is  encountered  in  attempting  to 

solve  the  deionization  equations  analytically  is  due  primarily  to  un- 
!  equal  rates  for  the  dissociative  recombination  and  mutual  neutrali¬ 

sation  reactions  —  if  these  reactions  are  assumed  to  proceed  at 
equal  rates,  the  deionization  equations  may  be  solved  exactly  for 
t  certain  source  functions.  However,  most  of  the  current  estimates 

of  the  rate  coefficients  give  values  for  these  two  reactions  which 
differ  by  factors  ranging  from  about  4  to  10.  Hence  it  is  necessary 
to  have  available  solutions  which  provide  for  this  difference. 

The  need  for  such  solutions  became  acute  in  preparing  the  revised 
edition  of  the  Electromagnetic  Blackout  Handbook  [Ref.  2].  The  large 
number  of  calculations  which  were  necessary  to  produce  graphical 
material  for  the  Handbook  made  the  use  of  numerical  methods  of 
solution  prohibitively  expensive.  On  the  other  hand,  the  approximate 
analytic  solutions  in  use  at  that  time  were  known  to  have  poor 
accuracy  in  certain  regimes.  It  was  apparent  that  a  need  existed 
for  new  approximations  giving  higher  accuracy  but  retaining  the 
essential  characteristic  of  rapid  computation. 

The  derivation  of  the  approximate  analytic  solutions  which  were 
used  in  preparation  of  material  for  the  Handbook  is  presented  in 
this  report.  In  addition  to  the  assumption  that  the  reaction  rate 
coefficients  are  constant  at  any  given  altitude,  there  are  tw» 


*  See  the  note  at  the  end  of  this  report  for  reference  to  a  recent 
Lincoln  Laboratory  report  which  describes  a  computer  solution 
employing  numerical  integration  methods. 
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approximations  inherent  in  the  derivation:  both  the  ionization 
source  strength  and  an  effective  positive  ion  recombination  coefficient, 
which  will  be  defined,  are  assumed  independent  of  time.  The  effect 
of  these  approximations  will  be  illustrated  an  will  methods  of 
reducing  the  inaccuracy  they  produce.  It  will  be  seen  that  although 
the  solutions  can  be  used  for  an  occasional  hand  calculation,  they 
are  best  suited  for  rapid  machine  computation.  A  Fortran -coded 
computer  program  employing  the  solutions  is  available  [Ref.  3]. 
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DERIVATION  OF  THE  SOLUTIONS 


The  differential  equations  which  describe  the  three- species  system 
may  be  written  as 


dN 

— ~  =  q  -  a  N  N  -  AN  +  DN 
dt  d  e  +  e 


dN 

— =  -a.N  N  +  AN  -  DN 
dt  1  -  +  e 


dN. 


— -  =  q  -  a.N  N,  -  Ot^N  N, 
dt  1  -  +  de  + 


where 


(1) 

(2) 

(3) 


N^  =  number  density  of  electrons 
N  =  number  density  of  negative  ions 
N+  =  number  density  of  positive  ions 

q  =  rate  of  ion  production  (source  strength) 

a  =  dissociative  recombination  reaction  rate  coefficient 
d 

a.  =  mutual  neutralization  reaction  rate  coefficient 
i 

A  =  attachment  rate  coefficient 
D  =  detachment  rate  coefficient 
These  equations  satisfy  the  charge -balance  relation 


N  =  N  +  N 
+  e 


(4) 


4 


RM  65TMP-7 


If  we  define  an  effective  positive  ion  recombination  coefficient  as 


a  = 


N  a.  +  N  a, 
-  i  e  d 


then  the  positive  ion  differential  equation  becomes 


(5) 


—  =  q  -  <*V  (6) 

We  now  introduce  the  approximations  that  q  and  QL  are  independent 
of  time,  in  which  case  the  solution  of  (6)  is  found  to  be 


l+Ce-2/^‘ 


1  -Ce-2/irt 


(7) 


where 


C  = 


N  (°)  + 

+  v  a 


(8) 


Using  (4)  to  eliminate  N_  from  (1)  and  substituting  (7)  for  N+,  the 
electron  density  differential  equation  becomes 


dN 
_ e 

dt 


1  +Ce 


-  2v/qat 


1  -  Ce 


ZVqpt  t 


-  A+D+a 


(9) 


This  is  a  linear  equation  with  variable  coefficients,  and  the  solution 
may  be  obtained  in  integral  form  through  classical  methods: 
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-  (A+D+  ~-S&y  (a+d+  ~-/<j5Qt 


N  (t)  = 


-2/qa' 


«d/a 


f- 


q(l-Ce-2^') 


(Hj/ot 


dt 


+  K 


(10) 


where  K  is  a  constant  to  be  determined  from  the  initial  conditions. 


The  integral  in  (10)  may  readily  be  obtained  in  closed  form  if  a  “  ad« 
The  result  is  sometimes  termed  the  "equal-alpha"  solution*  The 
equal-alpha  solution  applies  therefore  to  cases  in  which  the  dis¬ 
sociative  recombination  and  mutual  neutralization  reaction  rate 
coefficients  are  equal,  and  to  cases  in  which  the  negative  ion  popu¬ 
lation  is  negligible  with  respect  to  the  electron  (and  positive  ion) 
population.  This  latter  situation  arises  when  the  rate  of  electron 
attachment  is  small  relative  to  the  rate  at  which  negative  ions  are 
lost  through  detachment  and  mutual  neutralization.  Since  this  occurs 
in  the  daytime  ionosphere  above,  say,  80  km  or  at  roughly  the  same 
altitudes  in  the  nighttime  ionosphere  for  fairly  large  ion  production 
rates,  it  is  seen  that  the  equal-alpha  solution  can  be  essentially ‘an 
exact  solution  of  the  three-species  system(except  for  the  constant- 
q  approximation)  for  certain  cases  of  interest  even  when  a^  and  ad 
are  not  assumed  to  be  equal. 

When  a  does  not  coincide  with  a^,  the  integral  in  (10)  is  not  so 
easily  obtained.  Applying  the  binomial  expansion  to  the  troublesome 
factors  and  grouping  terms,  the  integral  in  (10)  may  be  put  in  the 
following  form: 


*  The  well-known  result  for 


Ne(t)  =V1  -Ce 


-  z/qa  t  ' 1 


-(A+Dv'tf7)  t 


Wfr 

A+D+yqa  l 


1  -  e 


-2/qOrt  ^  -(A  +  D+/qo$ 


(11) 


i 
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Je(A  +  D  +  k/qS)t  |q  +  Dy3_  .  (^kq  +  lk-ZlDyiJe-^4 
f  (k-U^qtlk-^D/^-je'4^' 


+  -£-  (k-l)(k-E)(k-3)[kq  +  (k-8)Dy|- 


/3L1  -8^qat 


-•••} 


dt 


where 
k  = 


(12) 


Of  I  2  VqQ»  t| 

This  expansion  io  convergent  for  |  Ce  1,  which  is  always 

true,  and  hence  may  be  integrated  term  by  term.  After  evaluating 

the  constant  K,  the  resulting  solution  for  Ng  is 


Nfi(t)  =(l  -Ce-2^®^  Kj(l  -C)kNe(0)e’(A+D+k^)t: 
q+P/o~  r.  _  -(A+D+k/^)t-i 

A+D+k/qoi  L  6  J 

kq+(k  Z)dJ~^  -2/qat  -(A+D+kyqaJt] 

- c  A+D+(k:2)^sr  Le  -e  J 

C3  ..  kq+^k  4^Dv/a^  r  -4/qat  -(A+D+k/qa)tl 

+  1T  Le  'e  J 

c3  (l  mi.  -i  kq+*k‘6,D/a  r  -6755 1  -(A+D+k/55)tl 

(k-D(k-E)  A+F+oTTiTS  Le  -e  J 


+  .  .  . 

+  IzCl! 

nf 


(k-1  )(k-2).  .  .  (k-n+1 ) 


kqt(k-2n  )Dft  .Zn/&t_  -(A+D+k7q5,f) 
A+D+(k-2n) /qa  |^e  e  J 


+  .  . 


(13) 
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It  may  be  seen  that  for  k  =  1,  equation  (13)  reduces  to  the  equal-alpha 
solution,  equation  ('ll),  as  would  be  expected. 

The  same  procedure  used  to  obtain  this  result  may  also  be  applied 
to  the  negative  ion  density  differential  equation  (2).  Usually,  how¬ 
ever,  an  adequate  solution  for  the  negative  ion  density  may  be  ob¬ 
tained  through  the  charge -balance  equation  (4).  This  can  result  in 
some  numerical  inaccuracy  when  the  positive  ion  and  electron  con¬ 
centrations  are  nearly  equal,  but  in  such  cases  there  is  usually 
little  need  for  great  accuracy  in  the  negative  ion  density  solution. 

Equations  (7)  and  (13)  are  general  solutions  to  the  three-species 
system,  subject  only  to  the  two  approximations  inherent  in  their 
development.  Before  turning  to  an  examination  of  the  effect  of  these 
approximations,  it  is  worthwhile  to  consider  some  potential  compu¬ 
tational  difficulties  which  may  arise  in  applying  these  results  and 
how  the  difficulties  can  be  avoided. 

In  computing  electron  density  from  (13),  enough  terms  are  included 
in  the  series  until  the  contribution  of  succeeding  terms  becomes 
negligible.  Usually  this  requires  only  a  few  terms,  although  the 
exact  number  needed  in  a  particular  case  depends  cm  the  values  of 
the  parameters.  For  cases  in  which  the  term  D  /-*—  is  large 
compared  to  q,  which  may  occur  for  low  source  strengths  in  the 
daytime,  it  is  possible  to  observe  a  false  null  in  the  series  which 
should  not  be  interpreted  as  indicating  that  sufficient  terms  have 
been  included.  This  occurs  when  k  is  nearly  equal  to  2n,  where 
n  =  1,  2,3,...,  and  results  from  the  small  magnitude  of  the  term 
having  (k  -  2n)D  in  the  numerator.  Even  though  this  term  may 
be  negligible,  the  next  term  which  involves  (k  -  2n-  2)D^~  will 
probably  be  significant.  Thus  the  term  following  the  first  term  which 
is  negligible  should  be  examined  before  the  decision  is  made  to 
terminate  the  series. 

A  more  serious  difficulty  arises  if  use  is  made  of  (7)  and  (13)  in 
cases  where  the  source  strength  q  is  very  small.  Both  (7)  end  (13) 
become  indeterminate  for  vanishing  q,  and  hence  are  useless  for 
computational  purposes  in  such  cases.  To  handle  these  cases  it  is 
necessary  to  develop  additional  forms  for  the  solutions.  The 
solution  for  the  positive  ion  density  when  q  is  of  negligible 
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importance  may  be  obtained  either  by  taking  the  limit  of  (7)  as 
q — >  0,  or  by  solving  (6)  after  eliminating  q.  The  solution  is  found 
to  be 


N+(0) 

N+(t)  =  1  +«N+(0)t 


(14) 


A  new  solution  for  the  electron  density  may  now  be  obtained  after 
substituting  (14)  for  N  in  the  differential  equation.  Equation  (1) 
becomes 


a  n  (°)  \ 

AtD*  i  +  OfN^(O)  t  )  Ne  11 51 

Note  that  q  has  been  left  in  the  electron  density  differential  equation. 
This  is  done  since  cases  arise  in  practice  wherein  the  source  strength 
is  an  important  term  in  the  electron  density  equation  even  though  it 
has  a  negligible  effect  on  the  positive  ion  density.  This  situation  is 
most  likely  to  occur  after  a  large  initial  level  of  ionization  in  the 
nighttime  at  fairly  low  altitudes  where  attachment  rates  are  large. 

The  inclusion  of  q  in  (15)  does  not  lead  to  any  computational  diffi¬ 
culty  for  the  case  where  q  =  0.  As  before,  the  solution  of  (15)  may 
be  obtained  in  integral  form  through  classical  techniques: 


dN  DN  (0) 

_ e  _  _ + 

dt  q  1  +aN+(0)t 


-(A+D)  t  f 

»i  e  1 

Tp  (A+D)t!a 

r  ittd 

1  l  +aN+(0)thr 

Ne  1  ~  r  ia<*  1 

l  +aN  (0)t  hjr-  l 

,Je  r 

L  +  J 

+DN+(0) 

* 

[l  +aN+(0)tj^r-"1J  dt 

+K 

(16) 

As  was  the  case  previously,  the  integration  is  readily  performed  in 
closed  form  if  a  =  yielding  the  equal-alpha  electron  density 
solution  for  small  q.  *  When  a  and  Ofj  are  unequal,  equation  (16) 


*  The  solution  of  ( 1 6)  for  a  =  it 

Ne(0)e'(A+D)t  N+(0)[D(A+D)-qo]  [^-(A+D)! 


Ne(t)  =  l+crN+(0)t 


A+D 


1  - 


(A+D)s 

-  (A+D)t 
1  +aN+(0)t 


j_e-(A+D)t  ■ 
,  1  +aN+(0)t  , 


(17) 
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may  be  put  into  a  form  suitable  for  computation  through  repeated 
integration  by  parts.  Two  such  expansions  are  possible,  and  each 
has  application.  For  cases  such  that 


1  +  aN+(0)t  <  2 


«N+  (0) 
A+D 


the  following  convergent  expansion  is  useful: 


N_(0)  e 


(A+D)t 


N  (t)  =  , 

e  [l  +  «N+(0)t]k 


+ 


-  (A+D)  t 


ka 


[l  +  oN+(0)tl 


taft  -  DIAtD)  If  .  ...  e„-^D)  1 

k'(k+l)  a8  N+(0,  \\  “  ♦«  [1+aN+(0),j 


~!k 


I 


IF  [*r«r]  -r'. 

T 


-  (A+D)t 


[l+oN+(0)tJ 


EaNJO)]  ([1+“\<°>*T  • 


(k+2)(k+3)  LaN+(0). 


(A+D)t  \ 


[l  +  OtN^(0)tJ 


+ 


_ izif. _ r. . a+d.-i 

(k+2)(k+3).  .  .  (k+n+l)LoN+(0)  J 


+ 


-  (A+D)t 
e _ 

[l  +  aN+(0)t]k 


(18) 


For  cases  such  that 
1  +  aN+(0)t£  2 


«N+(0) 

A+D 
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the  foregoing  expansion  converges  rather  slowly  and  hence  is  not  very 
suitable  for  computation.  In  such  cases  the  following  semiconvergent 
expansion  may  be  used  to  advantage: 


Ne(t)  * 


N  (0)  e 
_ e _ 


■(A+D)  t 


[l  +  aN+(0)tj 


A  +  D 


I 1-  -s 


-  (A+D)  t 


+  N  ,o,  DtA+SLJaa  J 

+  V°'  (A+D)a  ' 


[l  +aN+(0)tJk 

-{A+D)t 


l  +aN+(0)t 


ttN.(O) 


-  (k-1) 


+^.  _ 1 

A+D  |1+«N 


l  +«N+(0)t]1 


-(A+D)  t  \ 


[l  +  aN+(0)tJ3  [l  +  OtN+(0)  tj"  | 


+  (k-l)(k-2) 


aN+(0) 

a 

.  j  e-(A+D)t 

A+D 

[l  +  orN+(0)  tj  [l  +  aN+(0)  t]k 

„n 


+  (-  1)  (k-l)(k-2).  .  .  (k-n) 


aN+(0) 

n 

/  i 

A+D 

l  [i  +  aN  ,  (0)  t]n+1 

-(A+D)  t 


[l  +  aN+(0)t]] 


(19) 


This  series  should  be  terminated  when  either  or  both  of  the  factors 
(k-n)  or 


1 


-(A+D)  t 


^[i  +aN+(0)t  n+1  [l+«N+(0)t 


change  sign.  It  may  be  seen  that  for  k  =  1  equation  (19)  reduces  to 
the  equal -alpha  solution  valid  for  small  q,  equation  (17). 
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The  decision  to  use  equations  (14)  and  (18)  or  419)  in  place  of 
equations  (7)  and  (13)  may  be  made  by  first  obtaining  from  (14)  a 
lower  limit  of  N+,  i.  e.  ,  that  value  of  positive  ion  density  which 
would  result  if  there  were  no  continuing  source  of  ionization.  If  the 
ion  production  rate  term  in  the  positive  ion  differential  equation  (6) 
is  insignificant  compared  to  the  term  aN  2  using  the  lower  limit  of 
N+,  then  the  value  of  N+  obtained  from  (14)  is  the  proper  value  and 
either  (18)  or  (19)  should  be  employed  td  find  the  electron  density. 
Stated  mathematically,  equations  (14),  (MB),  arid  (19)  are  valid  when 


q 


«  a 


N+(0) 


1  +  aN+(0)t 


s 


otherwise  equations  (7)  and  (13)  should  be  used. 


One. other  computational  difficulty  can  occur  if  (13)  is  employed  for 
cases  involving  essentially  zero  initial  conditions,  low  values  of  q, 
and  relatively  small  values  of  time.  In  such  cases  the  initial  terms 
in  the  series  in  (13)  may  be  several  orders  of  magnitude  larger  than 
the  final  answer,  convergence  is  quite  slow,  and  numerical  in¬ 
accuracy  may  result.  For  this  type  of  problem  suitable  approxi¬ 
mations  may  be  made  to  arrive  at  a  closed-form  expression  for  the 
electron  density  solution.  Thus  for  the  conditions 


/cja  t  «  1 . 

we  may  make  the  following  approximations 

v/q^t  -/qoTt 
e  -  Ce  a»  2 


(20) 


/qat  _  -v/qot  _ 
e  H  +  Ce  %  2 


N+(0)  +  qt 


A 


Introducing  these  into  the  positive  ion  density  solution  (7)  there 
results 


N  (t)  *  N  (0)  +  qt 

T  T 


(21) 
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The  integral  form  of  the  electron  density  solution,  equation  (10), 
reduces  upon  substitution  of  the  above  approximations  to 


N  (t)  w  e 
e 


-(A+D)  t 


(A+D)  t 


(q  +  E>[N+(0)+  qt]^)dt  +  K 


(22) 


i  j 

Performing  the  integration  and  evaluating  the  constant  in  terms  of 

the  initial  condition  we  obtain 


N  (t)  N  (0)  e 
e  e 


-(A+D)t 


+ 


’  q  +  DN+(0) 

Dq  If,  -(A+D) t 

A+D 

-  (A+D)=>jre 

+  Qqt.,. 

A  +  D 


(23) 


Sufficient  relations  have  now  been  developed  to  handle  all  situations 
commonly  encountered  in  practice.  It  still  remains,  however,  to 
provide  some  procedure  for  determining  the  value  of  CL  to  be  used 
in  the  equations.  The  definition  of  a,  equation  (5),  is  not  directly 
usable  for  this  purpose  since  it  is  expressed  in  terms  of  the  solutions 
of  the  equations,  although  the  initial  value  of  a  can  be  found  from 
(5)  using  the  known  initial  conditions: 


0(0)  = 


a.N  (0)  +  a  ,N  (0) 
i  - _ d  e 

N+(0) 


(24) 


As  t - >»,  a  approaches  a  steady- state  value  which  may  be  found 

by  introducing  the  steady- state  ion  density  solutions  into  (5).  For 
this  purpose  it  is  convenient  to  rewrite  the  definition  of  a  in  terms 
of  the  ratio  of  negative  ion  density  to  electron  density: 


$  In  practice,  the  term  "quasi-equilibrium"  is  often  used  instead 
of  " steady- state"  since  generally  the  source  strength  q  is  not  inde¬ 
pendent  of  time  as  has  been  assumed  here.  However,  it  often 
happens  that  q  varies  slowly  enough  with  time  so  that  after  the 
initial  transients  die  out  the  ion  densities  approach  a  condition  where 
they  are  essentially  in  equilibrixxm  with  q  at  any  given  instant. 
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where 


a 


Xflf.  +  of 

i 

X  +  l 


N 


e 


Under  steady- state  conditions  we  have  from  (2)  and  (7) 


(25) 


(26) 


Substituting  (27)  into  (25)  there  results 


a(°°) 


Aa.  +Da,  +  of. a .  . 

i  d  i  d  V 

A+D+a.  . 

i  vaH 


(27) 


(28) 


Equation  (28)  yields  a  cubic  in  Of  which  may  be  solved  directly,  or 
an  iterative  procedure  may  be  used.  A  good  starting  point  for  the 
iteration  is  provided  by  replacing  Of  in  the  right-hand  member  by 
(Aor  +  DOf^)  /  (A+D).  This  is  the  limiting  form  of  0f(«)  for  the  case 
of  vanishing  q.  The  values  obtained  through  iteration  converge 
rapidly  to  the  solution  of  (28),  and  often  the  very  first  value  furnishes 
an  acceptable  approximation  to  the  solution. 


Once  the  initial  and  final  values  of  Of  have  been  obtained  from  (24) 
and  (28)  respectively,  a  reasonably  good  approximation  for  Of  at 
intermediate  times  has  been  found  to  be 


Of (t)  to  o<(0)  e 


-  (A+D)  t 


+  0f(®)  jji  -  e 


(A+D) 


'] 


(29) 
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ACCURACY  OF  THE  SOLUTIONS 


The  accuracy  of  the  solutions  developed  in  the  preceding  section  is 
dependent  upon  the  validity  of  the  "constant-q"  and  "constant -a" 
approximations  in  any  given  situation.  For  those  cases  in  which  q 
and  a  exhibit  little  variation  over  the  time  interval  of  interest,  the 
solutions  are  found  to  be  quite  accurate.  For  example,  if  and 
are  essentially  equal,  or  if  the  negative  ion  density  is  initially 
small  relative  to  the  electron  density  and  remains  so  due  to  low 
rates  of  attachment,  the  constant -a  approximation  is  valid  and  the 
accuracy  of  the  solutions  is  dependent  upon  the  validity  of  the 
constant-q  approximation  alone.  The  latter  situation  includes  many 
of  the  cases  encountered  in  the  analysis  of  nuclear  burst-produced 
ionization  at  altitudes  above  about  80  km. 


Another  type  of  problem  for  which  the  constant -a  approximation  is 
valid  is  one  in  which  sufficient  time  has  elapsed  to  attain  a  condition 
of  qua  si -equilibrium  with  a  slowly  varying  (or  constant)  source. 

In  this  case  a  stabilizes  at  the  value  given  by  equation  (28)  and  is 
essentially  constant.  The  quasi -equilibrium  term  in  equation  (13), 


N  (®) 
e 


A  +  D  +  a  , 
a 


(30) 


then  yields  a  result  which  is  an  exact  solution  for  the  quasi-equil¬ 
ibrium  case  providing  a  is  computed  exactly  from  equation  (28). 
Since  a  good  solution  to  (28)  is  readily  obtained  using  a  very  few 
iterations,  quite  acceptable  accuracy  (on  the  order  of  a  few  percent 
or  better)  is  realized  for  the  quasi- equilibrium  case  with  less 
computation  needed  than  by  most  other  methods  offering  comparable 
accuracy. 


A  determination  of  the  accuracy  of  the  solutions  for  problems  in¬ 
volving  altitudes  below  about  80  km  and  times  prior  to  the  establish¬ 
ment  of  quasi -equilibrium  is  somewhat  complicated  by  the  fact  that 
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a  number  of  parameters  are  involved,  each  of  which  may  range  over 
a  wide  spectrum  of  values.  The  accuracy  in  any  given  situation  is 
a  function  of  the  values  of  the  parameters.  To  evaluate  the  approxi¬ 
mate  solutions,  a  number  of  cases  have  been  computed  and  the  re¬ 
sults  compared  with  solutions  obtained  through  numerical  integration 
of  the  deionization  equations.  *  In  general,  the  maximum  error 
attributable  to  the  constant -a  approximation  occurs  during  the  attach¬ 
ment  phase  for  cases  in  which  the  initial  ion  concentrations  are  large, 
the  source  strength  is  relatively  small,  attachment  is  a  significant 
electron  loss  process,  and  the  dissociative  recombination  and  mutual 
neutralization  reaction  rate  coefficients  are  appreciably  different. 

The  attachment  phase  may  be  defined  as  that  period  of  time  when 
terms  involving  e"(A+D)  t;  acel important.  During  this  period  the 
true  value  of  a  experiences  a  pronounced  variation  with  time.  At 
later  times  a  is  a  slowly-varying  quantity  at  worst,  and  hence  the 
constant -a  approximation  becomes  quite  good. 

Consider,  for  example,  a  problem  defined  by  the  following  para¬ 
meters: 


N  (0)  =  N  (0)  =  1  x  10l(r  cm"3 
e  + 

q  =  1  x  104  cm*'5  sec"1 

h  =  50,  60,  70,  80  km  (day  and  night) 

The  reaction  rate  coefficients  for  these  altitudes  are  taken  from 
Reference  4  using  the  1962  U.S.  Standard  Atmosphere  for  atmos¬ 
pheric  density  and  temperature  values.  The  coefficient  values  are 
listed  in  Table  1 . 


*  The  numerical  integration  was  performed  using  a  computer 
program  written  by  C.  Luehr,  TEMPO.  This  program  uses  the 
SHARE  library  routine  R  W  INT  aft  its  basic  component.  This  is  a 
Fortran  version  of  R  W  -DE2F  which  integrates  a  system  of  N  simul¬ 
taneous,  first  order,  ordinary  differential  equations  with  the  option 
of  using  either  4th  order  Runge-Kutta  method,  or  4th  order  predictor- 
corrector  method  (Adams -Moulton)  with  either  fixed  or  variable 
step  size.  Double  precision  is  used  internally  to  control  round-off 
errors.  The  predictor-corrector  method  with  variable  step  size 
was  used  to  obtain  the  results  presented  in  this  report. 
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Table  1.  Reaction  rate  coefficients  used  in  calculations. 


h 

(km) 

(cm'J  sec" ') 

ot\ 

(cm^  sec-1) 

A 

(sec"1) 

D  (day) 
(sec"1) 

D  (niaht) 
(sec"1) 

50 

2,586x10'? 

6.477xl0’8 

5.842X101 

0.4415 

1.473xl0“3 

60 

2.738x  10’7 

6.384x  10"8 

4.954x  10° 

0.4408 

7.984X10"4 

70 

3. 188x  10'7 

6.795x  10'8 

3.577xl0-1 

0.4497 

1.801  xlO'3 

80 

3.875x  10"7 

7.454x  I0~8 

1 .687x  10"2 

0.4957 

2.050xl0"3 

The  electron  density  time  history  for  this  problem,  obtained  through 
numerical  integration,  is  presented  graphically  in  Figures  1  and  2 
for  daytime  and  nighttime  values  of  detachment,  respectively. 
Numerical  values  of  electron  density  at  times  of  0.  01,  0.  03,  0.  1, 
0.3,  1 . 0,  3.  0,  10,  30,  and  100  seconds  are  listed  in  Tables  2  and 
3.  The  first  column  for  each  altitude  contains  the  values  given  by 
numerical  integration.  These  are  considered  exact  for  the  purpose 
of  evaluating  the  approximate  analytic  solutions.  The  second  column 
for  each  altitude  contains  values  calculated  from  the  approximate 
solutions  derived  in  the  preceding  section. 

The  comparison  in  Tables  2  and  3  shows  that  the  approximate 
solutions  exhibit  relatively  little  error  at  any  time  at  the  80-km 
altitude,  a  conclusion  previously  deduced  from  the  fact  that  a  re¬ 
mains  essentially  constant  at  a  value  near  at  such  altitudes 
since  relatively  little  build-up  of  the  negative  ion  density  occurs. 
What  error  that  does  occur  at  80  km  is  largest  for  the  nighttime 
case,  which  is  not  surprising  since  larger  negative  ion  densities 
are  formed  for  the  lower  nighttime  detachment. 

At  the  50,  60,  and  70  km  altitudes,  the  error  in  the  approximate 
series  solutions  is  seen  to  be  a  pronounced  function  of  time,  be¬ 
coming  quite  large  at  certain  times  and  then  falling  rapidly  to 
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Table  2.  Comparison  of  results  for  N  (0)  =  N  (0)  =  10 10  cm  ,  q  =  104  cm 
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negligible  values.  The  time  of  maximum  error  is  a  function  of 
altitude,  and  by  referring  to  Figures  1  and  2  the  maximum  error  at 
each  altitude  is  observed  to  occur  during  the  period  when  the  electron 
density  is  decreasing  most  rapidly  with  time.  This  is  the  "attachment 
phase"  mentioned  previously,  and  the  error  is  a  result  of  a  significant 
variation  in  the  value  of  a  during  this  period.  Note  that  the  error  be¬ 
comes  quite  small  after  the  attachment  phase  has  ended,  and  by  100 
seconds  when  qua  si -equilibrium  has  essentially  become  established 
for  this  problem,  the  approximate  solutions  are  seen  to  be  very 
nearly  exact  at  all  altitudes. 

Since  the  approximation  which  is  responsible  for  the  error  consists 
of  the  assumption  that  a  is  a  constant  from  the  initial  time  to  the 
time  of  calculation,  we  might  expect  that  if  the  time  interval  over 
which  the  actual  value  of  a  exhibits  appreciable  variation  were 
divided  into  several  smaller  time  intervals,  and  the  results  at  the 
end  of  each  interval  were  used  as  initial  conditions  for  the  next, 
the  error  could  be  reduced.  Such  is  indeed  the  case,  and  the  amount 
by  which  the  error  is  reduced  is  often  substantial.  To  illustrate 
this,  sequential  calculations  were  performed  at  the  times  shown  in 
Tables  2  and  3,  beginning  with  one  earlier  time  of  0.  003  second. 

The  values  of  Ne  and  N+  obtained  at  each  time  were  used  as  Ne(0) 
and  N+  (0)  for  the  next  calculation.  The  results  obtained  by  this 
procedure  are  presented  in  the  third  column  for  each  altitude  in 
Tables  2  and  3.  In  every  case  the  maximum  error  is  lower  than 
for  the  separate  calculations,  often  significantly  lower.  Perhaps 
the  most  striking  comparison  occurs  for  the  60-km  nighttime  calcu¬ 
lation  at  0.  3  second.  The  separate  calculation  for  this  time  is 
over  2  orders  of  magnitude  low  (the  largest  discrepancy  found  for 
this  method  of  calculation),  whereas  the  value  obtained  by  calcu¬ 
lating  sequentially  is  only  5.  5  percent  low. 

The  f  oregoing  problem  was  chosen  for  illustration  since  the  combi¬ 
nation  of  parameters  it  involves  is  representative  of  something 
approaching  the  worst  case  which  is  likely  to  be  encountered  in 
practice,  as  far  as  the  accuracy  of  the  approximate  solutions  is 
concerned.  The  large  maximum  errors  found  in  this  example  are  a 
result  of  the  combination  of  the  following  two  conditions: 

1)  unequal  0^  and  0^  together  with  significant  build-up 
of  the  negative  ion  concentration  due  to  relatively 
large  attachment  rates,  causing  a  pronounced  variation 
of  a  with  time, 
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2)  large  initial  electron  and  positive  ion  concentrations 
together  with  a  relatively  low  source  strength,  causing 
the  recombination  terms  to  be  important  in  both  the 
positive  ion  and  electron  density  differential  equations 
and  hence  making  the  precise  value  and  time  history 
of  Ot  important. 

Tests  can  be  inferred  from  these  conditions  which  may  be  used  to 
predict  which  cases  are  likely  to  be  subject  to  significant  errors. 

Thus,  cases  which  satisfy  all  three  of  the  following  tests  are 
suspect: 

a  4  a(0) 

aNf(0)t»l,  (31) 

“rl 

(A+D)t  <  -jS- 

These  tests  may  be  applied  before  the  approximate  series  solutions 
are  employed  at  a  specific  time  in  a  particular  problem,  and  if  all 
three  are  satisfied  a  good  solution  may  be  obtained  by  computing  at 
one  or  two  earlier  times,  using  the  results  at  each  earlier  time  as 
initial  conditions  for  the  next  time  interval.  This  method  is  em¬ 
ployed  in  the  computer  program  described  in  Reference  3. 

While  the  maximum  inaccuracy  in  these  approximate  series  solutions 
may  appear  to  be  large,  and  indeed  it  is  for  single-time  calculations 
under  certain  conditions,  it  should  be  noted  that  over  much  of  the  time 
interval  of  interest  in  nuclear  weapons  effects  studies  the  error  is 
quite  small  in  this  problem,  and  the  solutions  will  be  found  to  be 
better  for  different  problems  involving  lower  initial  conditions. 

Another  point  to  be  considered  is  that  even  in  this  problem  the  series 
approximations  are  generally  much  better  than  other,  simpler  analytic 
approximations  which  are  commonly  used.  One  such  approximation 
is  the  equal-alpha  solution,  equations  (11)  and  (17).  The  equal-alpha 
solution  is  often  used  with  a  value  of  alpha  taken  as  (Afl^+Da^)/( A  +D). 
To  show  the  comparison  between  this  solution  and  the  series  approxi¬ 
mations,  values  of  electron  density  for  the  foregoing  problem  have 
been  obtained  from  the  equal-alpha  solution  using  the  above  expression 
for  alpha.  These  values  are  presented  in  the  fourth  column  for  each 
altitude  in  Tables  2  and  3.  Although  the  percent  error  comparison  is 
somewhat  misleading  since  the  equal-alpha  solutions  err  on  the  high 
side  while  the  series  solutions  exhibit  their  primary  inaccuracy  on  the 
low  side,  examination  of  Table  2  and  3  reveals  that,  in  general, 
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the  error  in  the  equal-alpha  solutions  is  greater  than  that  of  the 
series  solutions  and  persists  over  a  longer  time  period.  Further¬ 
more,  there  does  not  appear  to  be  a  simple  method  of  reducing  the 
error  in  the  equal-alpha  solutions,  whereas  the  error  in  the  series 
solutions  may  be  readily  reduced  by  the  sequential  calculation  pro¬ 
cedure. 

It  is  interesting  to  consider  a  second  problem  for  which  one  of  the 
tests  which  may  be  used  to  predict  regions  of  inaccuracy  in  the 
series  solutions  is  never  satisfied: 


N  (0)  =  N  (0)  0 

e  + 

q  =  1  x  104  cm"3  sec"1 

h  =  50,  60,  70,  80  km  (day  and  night) 

This  is  the  same  as  the  first  problem  considered  except  for  the 
initial  conditions.  The  comparison  of  the  numerical  integration 
results  with  those  obtained  from  separate  calculations  using  the 
series  solutions  is  shown  in  Tables  4  and  5  for  daytime  and  night¬ 
time  detachment,  respectively.  It  is  seen  that  for  this  type  of 
problem  the  series  solutions  are  good  to  within  a  few  percent  at 
all  times  and  all  altitudes.  No  need  exists  to  perform  sequential 
calculations  to  obtain  good  accuracy  in  this  case;  the  solution  at 
any  time  is  essentially  unchanged  if  calculations  are  performed 
sequentially. 

The  discussion  thus  far  has  considered  the  effect  of  the  constant -OL 
approximation  for  several  combinations  of  parameters.  A  similar 
discussion  applies  if  the  source  strength  q  exhibits  significant 
variation  over  the  problem  time  interval.  In  such  cases  it  is  possible 
to  reduce  the  inaccuracy  resulting  from  the  constant-q  approximation 
by  performing  a  series  of  calculations  during  the  period  of  signifi¬ 
cant  q-variation,  using  the  results  obtained  at  each  time  as  initial 
conditions  for  the  next  calculation.  For  many  problems  encountered 
in  nuclear  weapons  effects  studies,  however,  the  source  strength, 
does  not  vary  rapidly  enough  with  time  to  warrant  special  treatment. 
Exceptions  to  this  wherein  sequential  calculations  may  be  necessary 
because  of  q-variation  include  cases  where  beta-ray  ionization  is 
present  at  some  times  but  not  at  others,  and  cases  involving  a 
relatively  rapid  variation  in  the  gamma-ray  source  due  to  changing 
atmospheric  mass  between  the  source  and  the  point  of  interest. 
Another  situation  wherein  it  may  be  necessary  to  perform  sequential 
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calculations  occurs  if  fireball  radiance  causes  the  total  detachment 
rate  coefficient  D  to  vary  appreciably  with  time. 

The  approximate  solutions  presented  here  are  intended  primarily 
for  rapid  machine  computation.  While  they  can  be  used  for  hand 
calculation,  their  principal  advantage  lies  in  the  fact  that  they  can 
be  used  to  obtain  solutions  of  higher  accuracy  than  that  obtainable 
from  simpler  approximate  solutions,  and  with  a  much  lower  ex¬ 
penditure  in  computational  time  than  that  required  by  numerical 
integration  methods.  When  applied  to  machine  computation,  it  is 
a  relatively  simple  task  to  program  the  machine  to  sense  regions 
wherein  the  accuracy  may  be  inadequate  and  automatically  perform 
a  few  calculations  at  earlier  times  to  improve  the  accuracy.  In 
this  manner,  the  solutions  appear  to  be  capable  of  producing 
acceptably  accurate  results  for  all  types  of  problems. 


27 


RM  65TMP-7 


NOTE  ADDED  IN  PROOF 


During  the  proofing  of  this  report  for  publication,  the  author  was 
privileged  to  receive  a  copy  of  a  Lincoln  Laboratory  report  des¬ 
cribing  work  of  a  similar  nature  carried  out  jointly  by  Mt.  Auburn 
Associates  and  M.  I.  T.  Lincoln  Laboratory.  The  report  "A  Model 
and  a  Computer  Program  for  Atmospheric  Deionization  at  Early 
Times  Following  Nuclear  Explosions"  by  Frederick  L.  Ek 
[Project  Report  PA-77  (BMRS),  24  December  1964]  is  concerned 
with  basically  the  same  problem  as  that  considered  here  —  mfethods 
of  obtaining  acceptably  accurate  solutions  to  the  system  of  differential 
equations  which  describe  the  complex  atmospheric  deionization 
processes.  As  is  the  case  here,  the  treatment  is  concerned  primarily 
with  the  three-species  lumped-parameter  system. 

Whereas  the  approach  here  is  one  of  obtaining  approximate  analytic 
solutions  so  as  to  avoid  the  necessity  of  employing  numerical  inte¬ 
gration  techniques,  the  method  described  by  Ek  is  to  use  numerical 
integration  to  solve  the  time -varying  portion  of  the  problem, 
switching  to  steddy  -state  analytic  solutions  at  times  where  these  are 
sufficiently  close  to  the  numerical  solution.  The  numerical  inte¬ 
gration  method  is  the  same  as  that  used  here  to  check  the  accuracy 
of  the  series  solutions  (the  same  SHARE  numerical  integration  sub¬ 
routine  is  employed).  One  feature  of  the  computer  program  des¬ 
cribed  by  Ek  is  that  both  the  source  strength  and  the  detachment  rate 
coefficient  are  allowed  to  be  functions  of  time. 

It  still  appears  that  the  series  solutions  developed  here  can  be  used 
to  fill  the  void  hitherto  existing  between  simple  " equal-alpha" -type 
approximations  on  the  one  hand,  and  precise  numerical  integration 
on  the  other.  The  fact  that  regions  of  inaccuracy  may  be  detected 
a  priori  and  the  solutions  easily  improved  is  believed  to  be  a  very 
real  advantage.  This  combined  with  the  relatively  insignificant 
increase  in  computer  time  required  by  the  series  solutions  over 
simpler  and  less  accurate  approximate  solutions  appears  to  make 
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the  former  particularly  well  suited  for  application  in  computer 
programs  \such  as  TEMPO'S  WEPH  [Ref.  3  ]  and  perhaps 
Lincoln  Laboratory'.s  BOOM. 
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